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ABSTRACT 

This paper presents a librational solution for evolutions of parameters averaged over 
a synodic period in mean motion resonances in planar circular restricted three-body 
problem (PCR3BP) with non-gravitational effects taken into account. The librational 
solution is derived from a linearization of modified Lagrange’s planetary equations. 
The presented derivation respects properties of orbital evolutions in the mean motion 
resonances within the framework of the PCR3BP. All orbital evolutions in the PCR3BP 
with the non-gravitational effects can be described by four varying parameters. We 
used the semimajor axis, eccentricity, longitude of pericenter and resonant angular 
variable. The evolutions are found for all four parameters. The solution can be applied 
also in the case without the non-gravitational effects. We compared numerically and 
analytically obtained evolutions in the case when the non-gravitational effects are the 
Poynting-Robertson effect and the radial stellar wind. The librational solution is good 
approximation when the libration amplitude of the resonant angular variable is small. 

Subject headings: Interplanetary dust - Mean motion resonances - Celestial mechanics 
- Poynting-Robertson effect. Stellar wind 


1. Introduction 


Oscillations are often present when a physical system is near a stable state. The gravity of a 
star and a planet (or a planet and a satellite) that move according to a solution of the two body 
problem can perturb the motion of a body with negligible mass (restricted three-body problem). 
The dynamic of the body with negligible mass includes in this case also the so called mean motion 
resonances. In a mean motion resonance a ratio of orbital periods of the two minor bodies oscillates 
near a ratio of two natural numbers. 


The perturbed motion of the body captured in the mean motion resonance can be studied 
using a disturbing function. The disturbing function is often expanded using Fourier series of the 
Laplacian type (e.g., Murray &: Dermott| 1999). Usability of this access is limited by two factors. 
First, a finite order of the expansion makes it practically unusable for large eccentricities. The 
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second limit is the approximation of the averaged disturbing function by a single periodic resonant 
term. Murray & Dermott (1999) substituted the truncated averaged Fourier series in the time 


derivatives of orbital elements given by Lagrange’s planetary equations. By integrating of these 
equations individually they obtained dependencies of orbital elements on time. To derive the time 
dependencies they assumed that the only time-varying quantities in the equations for the time 
derivatives of orbital elements are in the trigonometric arguments of the single periodic resonant 
term and that the longitude of pericenter increases linearly with time at a constant rate determined 
by secular theory. The obtained results were applied for an asteroid in the 1/2 interior resonance 
with Jupiter (Murray Sz Dermott||1999 p. 259). Another example of the time dependence obtained 


using the Fourier series expansion of the disturbing function can be found in Greenberg (1973). 


Lagrange’s planetary equations for the planar circular restricted three-body problem (PCR3BP) 
including tidal dissipation were solved simultaneously using several crude approximations. The 
solution of model was applied to the 4/3 exterior resonance between Saturn’s satellites Titan and 
Hyperion. It has not been possible to match exactly the behavior of Titan-Hyperion resonance using 
this model because of the approximations used ( Greenberg 1973). Even so, the model provided 
insight into the underlying mechanisms. 

For the mean motion resonances in the PCR3BP canonical equations can be written using 
Hamiltonian formalism and averaged ever a synodic period. This access gives a system with one 
degree of freedom that is integrable. However, the obtained system is difficult to solve and analytical 
solution in general case was not found. The dynamics of dust particles captured in the mean motion 
resonances and simultaneously affected by non-gravitational effects was frequently investigated 
during the last three decade^ Influence of electromagnetic radiation of a star on the dynamics of 


dust particles can be described by the Poynting-Roberson (PR) effect (Poynting 1904; Robertson 


1937; Klacka 2004; Klacka et al. 2014). In Beauge Sz Ferraz-Mello (1994) the equations of motion 


of a dust particle captured in a mean motion resonance in the PCR3BP with the PR effect were 


written in a near canonical form. Beauge Sz Ferraz-Mello (1994) transformed the near canonical 


equations to a system of equations suitable for searching of librational points. In order to determine 
the stability of a librational point they mentioned a linearization performed around the libration 
point and a solution of obtained characteristic equation of the system. The stability of the libration 
points was determined with a similar procedure also in Sidlichovsky &: Nesvorny (1994). In both 


previously mentioned papers specihc procedure used for solving of the linearzied system cannot 


be found. Even more, the linearized system in Beauge Sz Eerraz-Mello (1994) and Sidlichovsky & 


m 


Nesvorny (1994) is different. In Beauge Sz Eerraz-Mello (1994) the system has four equations and 


Sidlichovsky & Nesvorny (1994) has only three equations. 


In this paper we derive the solution of the linearized system of equations describing evolutions 


^For example we can mention the following papers: 

Jackson & Zook 

(19891 

Weidenschilling & Jackson 

(I993I); 

Marzari & Vanzani (19941; Lion & Zook (19951; Lion 

et al. 

(|1995l; |Liou & Zook| (|1997 |1999l; |Moro-Martm & 

Malhotral (|2002|); |Holmes et al.| (|2003|); |Kuchner & Holman ( 

2003|; Deller & Maddison| (|2005|); Stark & Kuchner 

(20081; Reach (20101; Mustill & Wyatt 

(20111; Pastor (20131; 

Lhotka & Celletti (20151; Shannon et al. (20151. 




































































































































































of dust particles captured in mean motion resonances in the PCR3BP with non-gravitational effect. 
We apply the derived solution in the case when the non-gravitational effects are the electromagnetic 
radiation of a star and its radial stellar wind. 


2. Averaged resonant equations 


A secular evolution of a dust particle in an orbit around the Sun is significantly influenced also 
by non-gravitational effects. In order to describe the secular evolution in the most of practical cases 
time derivatives of orbital elements averaged over some time interval can be used. The averaged 
time derivatives of orbital elements of the dust particle’s orbit caused by the non-gravitational 


effects can be calculated using the Gaussian perturbation equations of celestial mechanics (Danby 


1988; Murray Sz Dermott 1999). The secular time derivatives determine the secular evolution in 


the case without the gravitational influence of a planet. In order to obtain a system of equations 
describing the secular evolution in a mean motion resonance with a planet we use Lagrange’s 


planetary equations (Brouwer &: Clemence 1961 Danby 1988). The non-gravitational effects for 


which is possible to calculate the secular time derivatives of orbital elements can be included in 
Lagrange’s planetary equations. In a planar case (when the dust particle’s orbit lies in the planet 
orbital plane) modified Lagrange’s planetary equations which include the non-gravitational effects 


are 


dcTh 

dt 


da 

dt 

de 

dt 

dCj 

dt 

dn 


+ 




da \ 

dt / gp 

a BRq 
Le duj 
doj \ 

dt / pp 
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+ 


de \ 
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+ 


/ doh dn\ 

+ \ dt + dt/pp 


( 1 ) 


Here, a is the semimajor axis of the particle orbit, e is the eccentricity, u is the longitude of 
pericenter, and the angle fib is defined so that the mean anomaly can be computed from M = nt 
+ (Tb ( Bate, Mueller fc White||l971 ). Where n = V^/r(l - 13)/a^ is the mean motion of the particle, 
/i = GqM^, Gq is the gravitational constant, and M* is the mass of the star. The parameter /3 
is defined as the ratio between the electromagnetic radiation pressure force and the gravitational 
force between the star and the particle at rest with respect to the star. L = ^yJL{^^^~P)a. Rq is 
the disturbing function 

1 r • rp 


Rg = GoMp 


( 2 ) 


|r — rp| Tp 

In the disturbing function: Mp is the mass of the planet, rp is the position vector of the planet 
with respect to the star, rp = |rp|, and r is the position vector of the dust particle with respect to 
the star. The partial derivative of the disturbing function with respect to the semimajor axis in the 
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last equation in Eqs Q is calculated with an assumption that the mean motion of the particle n 
is not a function of the semimajor axis. All terms in Eqs. Q represent quantities averaged over a 
synodic period. {da/dt)EF, {de/dt)EF, {du}/dt)EF, and (dcib/dt + t dn/dt)EF are caused by the non- 
gravitational effects only. Eqs. 0 already describe the secular evolution of the dust particle which 
is captured in the mean motion resonance and simultaneously affected by the non-gravitational 
effects. For the study of a specific mean motion resonance, it is convenient to define a resonant 
angular variable (Beauge Sz Ferraz-Mello 1993, 1994; [Sidlichovsky &: Nesvorny 1994; Pastor 2014) 


p + q. , ~ 

cr = -Ap — s\ — u , 


( 3 ) 


where p and q are two integers (resonant numbers), Ap is the mean longitude of the planet, A = M 
+ d; = nt + Ub + w is the mean longitude of the dust particle, and s = pjq- In the mean motion 
resonance a is librating rather than circulating. The resonant angular variable is not yet averaged. 
For the time derivative of a we obtain from its definition 


da 

dt 


p + q 

-np 

q 


— sn — s 


/ dah dn duj\ 
\ dt dt dt ) 


dw 

dt 


( 4 ) 


The above equation can be averaged over the synodic period. If we use in the averaged result the 
last two averaged equations in Eqs. Q, then we get 


da 

dt 


" n I n M , 2sL 

p + q / dui\ / d(Tb ^ \ 

/ EF \ dt dt j gp 


dRo p^q 

+ np-ns 


( 5 ) 


The averaged disturbing function Rq^ can be rewritten as a function of the variables a, e, a), and 
a. In this case the following relations hold 

dRo _ dRo, 
da\i ^ da ’ 

ORg _ p + q ORg . . 

did ~ q da ' ^ ^ 

We can use in the first two equations in the system of equations given by Eqs. 0 relations Eqs. 
0. The last equation in Eqs. 0 can be replaced with equivalent Eq. 0. By this we obtain 
a system that enables to study the secular orbital evolution of the dust particle captured in the 
specific mean motion resonance given by the resonant numbers p and q under the action of the 
non-gravitational effects, 
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( 7 ) 


The last term in the angle brackets despite of its complicated meaning can be straightforwardly 
obtained using ( Bate et al.||1971 ) 

o 2 r 


/ + t—\ - / — 

\ dt dt j j^p 


na 


Or 


cos/ 


1 + e cos / 


ap - 


sin / 2 + e cos / 
e 1 + e cos / 


( 8 ) 


where / is the true anomaly, ap and op are the radial and transversal components of the acceleration 
caused by the non-gravitational effects. 


The system of equa tions given by Eqs. ([7|) is different from systems considered in Beauge & 
Ferraz-Mello (1994) and Sidlichovsky & Nesvorny (1994). The system in Beauge & Ferraz-Mello 


(1994) does not take into account the non-gravitational effects for which the secular variations of 


the particle orbit depend on the orientation of the orbit in space. The system in Sidlichovsky & 


Nesvorny (1994) contains only three equations. 


3. Linearization of averaged resonant equations 

No general method exists for solving nonlinear differential equations in the system Eqs. 0. In 
practice, the best that can be accomplished is to study a linearization based upon initial conditions 
for the function and its derivatives. In the vicinity of an initial point oq, cq, wq, and a^ we use 
notation 


5a — a Oq , 

5e = e — eo , 

5ui — O! Coq , 

5a ^ ^0 • 

The time will be measured from an initial time to- Hence 

5t = t. 


(9) 


( 10 ) 


The linearization of averaged resonant equations (Eqs. finally gives 
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d5e 
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The terms in parenthesis before the deltas on the right-hand sides i 
we can simply white 


m 


( 11 ) 

Eqs. © are constant therefore 


6a = Aba + Bbe + Cbcj + Dba + Et -|- F , 

be = Gba + H be + I bijj + Jba + Kt + L , 

b^ = Mba + Nbe -t- Ob^ + Pba + Qt -|- R , 

ba = Sba + Tbe + Ub^ + Vba+Wt + X . ( 12 ) 


This system describes solution of system Eqs. Q during a short time interval after the initial time. 
It is possible to obtain an equation for one chosen delta by an elimination of the remaining deltas 
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using all equations in the system. The obtained equations for the separated deltas are 


5 a + Aa 3 (5 a + ^a2 '^a + AqI <5^ + A^q + ^at ^ + A^ = 0 , 

5 e + Ae 3 (i e + ^e2 + Agl jg + AgO <5g + Agj t + Ag = 0 , 

5 d; + A^3 5 G) A(I,2 + A(2,i 5cj + A^q ^Cj + ^u)t t + A^ = 0 , 

(j o- + A(j 3 j (t + A(j 2 5fT + Ag-i (5(j + Ao-0 <^0- + Ao-i t + Ag- = 0 . (13) 


One possible way how we can obtain the separated equation for 6a in Eqs. (13) is to calculate the 
following time derivatives of the first equation in Eqs. (12). 


6 a = A 6 a + B 6 e +0(5(2, + D 6 a + Et + F , 

6 a = OL 26 a + /52(5g + 72<5d; + 626 a + 62 ^ + C 2 , 

6 a = 0 : 36 a + /53<5g + 73<5d; + + ^3^ + C 3 , 

6 a = Oi46a + /34(5g + 74^ + ^ 46 ^ + £ 4 ^ + C4 , (1^) 


here ai, ^i, ji, 6i, e;, and Q for / = 2, 3, 4 are determined by the constants in Eqs. 


We can 


substitute Eqs. (14) in the first equation in Eqs. (13). Now, when we realise that the first equation 
in Eqs. (13) should be valid for arbitrary deltas, we obtain the following system of equations 


0(4 + Aa3 0:3 + Aa2 02 + AqI A + A^O — 0 , 

Pa + A(j3 P3 + Ka2 P2 + Aal i? = 0 , 

74 + Aa3 73 + Aa2 72 + A^l C = 0 , 

<54 + Aa3 63 + Ka2 62 + A^i H = 0 , 
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C4 + Aa3 C3 + Aa2 C2 + A^l F + A^ = 0 . (15) 


The solution of Eqs. (15) gives unknown Deltas in the separated equation for 6a- In the used 
notation we obtain for all separated equations in Eqs. (13) 
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( 21 ) 


We can simplify the equations above using properties of the mean motion resonances in the 
PCR3BP. Each of the equations in Eqs. 0 can be written as a sum of two terms 


dt 


-^1(^0) + P2(w) 


( 22 ) 
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where 'h denotes one of the variables (a, e, w, a), Fi is produced by the gravitational forces in 
the considered system and F 2 is caused by the non-gravitational effects. For the non-gravitational 
effects with a rotational symmetry around the star ^ 2 ( 0 ;) is a constant. For partial derivatives of 
Eq. (22) with respect to a) and a we obtain 


From Eq. ^ we have 


da 

doj 


d 

dT 

dFi 

da 

dF2 


dib 

dt 

da 

dob 

dbb ’ 


d 

dT 

dF^ 

dFo dCb 

(23) 


dt 

da 

+ ■ 

dCb da 


P + ' 

q 

dCb 

q 

(24) 


Q 


da 

p + q 


Substitution of Eqs. (24) in Eqs. (23) leads to the relation 

5 d'h p + q d d'^ 


doj dt 


q da dt 


(25) 


This relation implies 


C= -P±^D, 

Q 

I = J, 

Q 

0 = -^P, 

Q 

U=-P±^V. 
d 


(26) 


Therefore, determinants for Aq, and Ag* have linear proportionality between columns. In this 

(27) 


case properties of determinants reduce corresponding equations in Eqs. (19)-(20) to 

Aq = Aat = Let = 0 . 


General solution of Eqs. (13) with substituted Aq = 0 has form 

Ai A 2 As 2Ai A| 


(28) 


Here, are complex constants, B-qt are real constant numbers (as we will see later), and Aj, i = 
1, 2, 3 are all roots of a cubic equation with real coefficients 


A^ + AsA^ + A 2 A + Ai — 0 


(29) 


The roots of any cubic equation with real coefficients are always three real numbers or one real 
number and two complex numbers which are complex conjugate to each other. The roots of Eq. 
(29) can be written in the following form 

A = y , 


(30) 
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where u and v are the roots of quadratic equation 


y^ + 



A3A2 

3 



1 




(31) 


Eqs. (27) and (28) imply that the semimajor axis {6a) and the eccentricity {6e) cannot have terms 
varying quadratically with the time for this linearized system. The condition that the particle 
is in a mean motion resonance implies that the semimajor axis cannot have also a term linearly 
proportional to the time. Hence 



B 
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E 

Aa = - 

H 
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- 
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N 

O 
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T 

V 

W 


(32) 


where also Eqs. (26) were used. Any dependence between rows of the determinants in Eq. (32) is 
physically worthless. The dependence between columns in the determinants suggests that Eq. (32) 
is satisfied with identities 


E =^iB+ C 2 C + CsD , 
K =^,H + 6/ + 6^ , 

Q =^iN + ^20 + ^3P , 
W = ^,T + i2U + , 


(33) 


where are real numbers. Adding of the column proportional to the partial derivatives with respect 


to semimajor axis {A, G, M, S) on the right-hand side of Eqs. (33) would give Aq / 0 in Eq. (32) 
In this case the semimajor axis would have the term varying linearly with the time. Hence, Eqs 


(33) without these terms are consistent with the resonant behavior of the evolution of the semimajor 


axis. Eqs (26) and (33) cause that also the longitude of pericenter and the resonant angular variable 
cannot have terms varying quadratically with time (Eqs. [20] ). We come to conclusion that for the 
sought for solution of the linearized system must hold 


Aa* — Ap/ — Arjt — Aat — 0 


(34) 


Remaining equations in Eqs. (21) are simplified with substitutions of Eqs. (26) and (33) into the 
following form 
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= 


A 

G 

M 

S 


B 

H 

N 

T 


C 

I 

O 

U 


F 

L 

R 

X 


- 6 


A 

G 

S 


B 

H 

T 


C 

I 

U 


- ?3 
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G 

S 


B 

H 

T 


D 

J 

V 


-6 


A 

M 

S 


G 

O 

U 


B 

N 

T 


(37) 


Our solution of the linearized system allows the eccentricity, the longitude of pericenter and the 
resonant angular variable have the term varying linearly with time. However, after substitution 


of the solution (Eq. 28) in Eqs. (12) all terms varying linearly with time on the right-hand sides 


must cancel each other otherwise we would obtain on the left-hand sides of Eqs. (12) terms varying 
linearly with time. Thus, we can write 

E =^B + ^G + ^D , 

Ai Ai Ai 

K + ^J , 

Ai Ai Ai 

Q = ^N+^O + ^P , 

Ai Ai Ai 

W = ^T +^U + ^V . 


Ai 


Ai 


Eqs. (33) and Eqs. (38) will be satisfied when 


Ag A^^ 

If we use Eqs. (26) in Eq. (18), then we obtain for Ai 


Ai 


P — 


Ai = - 


A 

B 

C 


A 

B 

D 

G 

H 

I 

— 

G 

H 

J 

M 

N 

O 


S 

T 

V 


(38) 


(39) 


(40) 


A comparison of Eq. (35) with Eq. (40) yields that the first equation in Eqs. (39) is trivially 


satisfied for all ^i. The second equation in Eqs. (39) gives 

( 






+ 6 


A D B 
M P N 
S V T 


A B G F 
G H I L 
M N O R 
S T U X 

\ 


- ?2 


A B G 
G H I 
M N O 


- 


A B D 
G H J 
MNP 


(41) 


where also Eqs. (|26|) were used. The last equation is equivalent with 

A 


B G F 
G H I L 
M N O R 
S T U X 


-6 


A G B 
MON 
S U T 


- 6 


A B G 
G H I 
S T U 


+ 6 


A 
G 

M N O 


B G 
H I 


= 0 


(42) 
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If we use Eq. (37) in Eq. (42), then we obtain 



A 

B 

G 


A 

B 

D 

Ao- = -^3 

G 

H 

I 

- 

G 

77 

J 


M 

N 

0 


S 

T 

V 


= ^3^1 


(43) 


In other words, if ^2 = then ^3 = A^/Ai. Hence, the obtained solution of the linearized 


system is self consistent (see Eqs. 39). We have obtained three free parameters ^ 1 , ^2 and ^ 2 - These 


parameters determine unknown partial derivatives with respect to time in Eqs. ( 0 . Substituting 


Eqs. (38) in Eqs. (12) one obtains 


da 

dt 


A« 


Ai 

A, 


A, 


= Aa +B{e + ^t] + C { Co + ] +D[a + ^t] + F* , 


Ai 

A,: 


A„ 


^ = G a + 77 ^e + J + 7 ) + J ( cr + '-^t ) + L* , 

^=A7 a + iV^e+ + O ) + P ( cr + ] + R* , 

^=Sa +rfe+^t^ +u(co + ^t] + E ( cr + ) + X* . 


Ai 

Ai 

Ai 

A, 


dt \ Ai J \ Ai 

This suggests validity of the following equations 

^ /ci^\ d^/'Ae\ d_ \ d_ /c«r\ f A^ 

dt \ dt J de \ dt ) dt \Ai ) doj \ dt ) dt \Ai ) ^ da \ dt ) dt \Ai 


(44) 


(45) 


Said another way, the direct time dependence always appears at some of the variables e, Co and 
cr. We can speculate if this property belongs only to our linearized system or it is a heritage of 
full solution of the system given by Eqs. (HI- For completeness, Eqs. (|45|) are equivalent with 


Eqs. (38). The partial derivatives with respect to time in Eqs. 0 are unknown and cannot be 


determined (see also Eqs. 21). We can guess values to be negatively taken, over a libration period 


averaged, slope of the evolution of the eccentricity, the longitude of pericenter and the resonant 
angular variable. If this assumption is usable, we obtain oscillations of a variable modulated on a 
line that has the slope equal to the time derivative of the variable averaged over the libration period. 
If there is rapid advance in the longitude of pericenter during a capture of the dust particle in a 
mean motion resonance in the PCR3BP with the non-gravitational effects, then the most natural 
choice seems 

Ap „ Acj _ A„ 


w = 0 
Ai 


— dt , = 0 

dt ’ Ai 


(46) 


We have adjusted the properties of the sought for solution according to the properties of 
the evolutions of dust particles in the mean motion resonances with the influence from the non- 
gravitational effects taken into account. After omitting the zero terms remaining terms in the 


solved system are (compare to Eq. 13) 


5 a + A 3 (i a + A 2 5a + Ai (5a — 0 
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5 e + ^3 <5 e + ^2 + -^e = 0 , 

5 + A3 5 ii, + A2 5(2, + Ai (5i, + A,2, = 0 , 

5 (7 + A3 5 o- + A2 5o- + Ai (5o- + Ag- = 0 . 


The sought for solution of the system given by Eqs. (47) has the form 


5a = 4 ^e^i' + ^ ^ 

Ai A2 A3 

A — I x^t I ^e3 A3t _ I D 

+A 2 " +A 3 " A/+^^’ 

Ai A2 A3 Ai 

. 

Ai A2 A3 Ai 


(47) 


(48) 


In the next step we determine the complex constants Aai, Aei, A( 2 ,i and Ag-i as well as Ba, B^, B^ 
and Bfj from initial conditions. We have assumed that in time zero 5a(0) = 5e(0) = 5,^(0) = 5o-(0) 
= 0. Therefore 


A {r\\ Aal I Aa2 Aa3 . p 

^a(O) = -h -h -h .Ba = 0 , 

Al A 2 A3 

_ , , Api Apo Ap3 „ 

4(0) = ^+ ^+ ^+5e=0, 

Al A2 A3 

5(2(0) = + ^ + + 5(2 = 0 , 

Al A2 A3 

<5^(0) = 4^+ 4^+ 4^+5,, = 0. 

Al A 2 A3 


(49) 


Equalling of Eqs. (12) and the differentiation of Eqs. (48) with respect to time in the time zero 
gives 


5 a( 0 ) — Aal + Aa 2 + Aa 3 — F , 

5 e( 0 ) = Agi + Ae2 + Ae 3 — — = L , 

Al 

4 ( 0 ) = A (21 + Acj2 + A (23 — = R , 


5o-(0) — Aa-l + Act2 + A(j3 — — — X 

Al 


(50) 


Equalling of the differentiation of Eqs. (12) with respect to time and the second differentiation of 
Eqs. (48) with respect to time in the time zero gives 


5 a( 0 ) — AaiAi + Aa 2 A 2 + AasXs — AF + BL + CR + DX + E , 
5 e( 0 ) = TleiAi + 7 le 2 A 2 + Ae 3 A 3 = GF + HL + I R + JX + K , 
4 ( 0 ) = A( 2 iAi + A( 22 A 2 + A, 23.53 = MF + NL + OR PX Q , 
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5,(0) = ^,iAi + = SF +TL + UR+ VX + W . 


(51) 


5 e(0) — AelXf + ^62^2 + ^63^3 — 


Finally, equalling of the second differentiation of Eqs. (12) with respect to time and the third 
differentiation of Eqs. (48) with respect to time in the time zero yields 

5 a(0) = AaiXf + ^,2^2 + = A [AF + BL + CR + DX F E ) 

+ B {GF +HL + IR +JX +K) 

+ C {MF + NL + OR + PX +Q) 

+ D{SF +TL +UR + VX + W) , 

G {AF + BL +GR + DX + E) 

+ H {GF +HL + IR +JX +K) 

+ I {MF + NL + OR + PX +Q) 

+ J{SF +TL +UR + VX + W) , 

M{AF + BL +GR + DX + E) 

+ N {GF +HL + IR +JX +K) 

+ O {MF + NL + OR + PX +Q) 

+ P {SF +TL +UR + VX + W) , 

S {AF +BL +GR + DX + E) 

+ T {GF +HL + IR +JX +K) 

+ U {MF + NL + OR + PX +Q) 

+ V {SF +TL +UR + VX + W) . (52) 


6u){0) — AcjiXi + ^^> 2^2 + Aq^X^ — 


S (t( 0 ) — AaiXi + ^, 2^2 + ^0-3^3 — 


From equation Eqs. (50)-(52) Aai, Aei, Acji and A^i can be obtained. The obtained equations are 

5 5 (( 0 ) — 5iif(0) (A 2 + A 3 ) + 


= 


Am = 


^'I'3 = 


I 

+ Ai 


A 2 A 3 


(Ai — A2) (Ai - 

5 5,(0) — 5 <i,( 0 ) (Ai + A3) + I 

- A3) 

(•suo) + 

1 A1A3 

(Ai — A2) (A3 - 

5 '['(O) — 5 vi/( 0 ) (Ai + A2) + I 

- A2) 

(■i-W + If) 

1 A1A2 


(Ai — A 3 ) (A 2 — A 3 ) 


(53) 


Ai^i in Eqs. (53) are related to Xi in such a way that if Ai and A 2 are complex conjugate to each 
other and A 3 is a real number, then also Aij,i and A-q /2 are complex conjugate to each other and 
Aifs is a real number. This property holds for any permutation of not equal indexes i. From Eqs. 

gives 

(A 1 A 2 + A 1 A 3 + A 2 A 3 ) 

= -- — - • (54) 


(49) we can obtain Ba, Be, Bc^ and i?,. Eqs. 

5ii-(0) — 5 ^ 1 ( 0 ) (Ai + A2 + A3) + 


A 1 A 2 A 3 


Bi^ are always real numbers. 
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Stellar radiation and its influence on orbital evolution of cosmic dust 


In astrophysical modelling of the orbital evolution of dust grains it is often assumed that 
a spherical particle can be used as an approximation to real, the arbitrarily shaped, particle. 
Influence of electromagnetic radiation on the motion of a homogeneous spherical dust particle can 


be described using the Poynting-Robertson (PR) effect (Poynting 1904; Robertson 1937 Klacka 


2004; Klacka et al.|2014 ). The acceleration of the dust particle caused by the PR effect in a reference 
frame associated with the source of radiation (star) is 


— = rR- 
dt 


U • CR \ ^ V 

1 -fiR- 

c / c 


(55) 


where r is the radial distance between the star and the dust particle, cr is the unit vector directed 
from the star to the particle, v the velocity of the particle with respect to the star and c is the 
speed of light in vacuum. For the parameter (3 we have from its definition after Eqs. 0 


/3 = 


3L^Q' 


pr 


l67rcfj,Rdg 


(56) 


where L* is the stellar luminosity, is the dimensionless efficiency factor for the radiation pressure 
averaged over the stellar spectrum and calculated for the radial direction (Qpr = 1 for a perfectly 
absorbing sphere), and is the radius of the dust particle with mass density g. Expanding solar 
corona is continuous flux of solar wind inside the heliosphere. Stellar winds were also observed in 


the vicinity of other solar-like stars (Wood et al. 2002, 2005). The motion of the dust particles 


orbiting such a star can be affected by the stellar wind. It is possible to derive the acceleration 


caused by wind corpuscules impinging on the dust particle using the relativity theory (Klacka & 


Saniga 1993; Klacka et al. 2012). A radial stellar wind produces the following acceleration of the 


dust particle in accuracy to first order in u/c (u is the speed of the dust particle with respect to 
the star), first order in u/c {u is the speed of the stellar wind with respect to the star) and first 
order in v/u 

V ■ e^\ ^ V 


dv T] u g, 
dt c 


1 - 


u 


cr 


u 


(57) 


r] is to the given accuracy the ratio of the stellar wind energy to the stellar electromagnetic radiation 
energy, both radiated per unit time 


7 ] = 


dvrr^tt 




/ ^ ^SW J^SW jC 


(58) 


where rusw j and Usw j, 3 = ^ N, are the masses and concentrations of the stellar wind particles 


at a distance r from the star (u = 450 km/s and rj = 0.38 for the Sun, Klacka et al.||20f2 ). We are 
interested in the motion of a dust particle in the frame of reference associated with the star. When 
we add the gravitational accelerations from the star and the planet, we obtain the final equation of 
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motion of the dust grain in the PCR3BP with the electromagnetic radiation and the stellar wind 
in the reference frame associated with the star 


dv fj, GoMp 

^ = —2 “ P) ^ 

at \r — rpr 


(f- fp) - 


GqMp 







v-ep^ V 

-eR + - 

c c 


(59) 


At summation of Eq. (55) and Eq. (57) in Eq. (59) is assumed that {rj/Q'^^){u/c) <C 1. 


Using the last term in Eq. (59) as a perturbation of the orbital motion in a Keplerian potential 


—^(1 —/3)/r, we obtain from the Gaussian perturbation equations of celestial mechanics (e.g. Danby 
1988) the following averaged values in Eqs. ([^ 


da \ 
dt / 


EF 


7 ?- (2 + 3 e^) , 


de \ 

dt / Ep 
did \ 

dt / 



= 0 


EF 


/ _L +^\ 

\ dt ^ dt / 


= 0 


(60) 


EF 


The substitution of the partial derivatives of Eqs. (60) with respect to a, e, Cj and a in Eqs. (0 
determines corresponding coefficients in Eqs. (12). The coefficients are shown in Appendix [A} 


Numerical results 


In this section we compare analytical and numerical results. Real resonant libration for a 
spherical dust particle can be obtained from numerical solution of Eq. ( |59[ ) in any specific mean 
motion resonance. Equations in the system of equations given by Eqs. ([^ are averaged over the 


synodic period. All parameters substituted in Eq. (Al) should be averaged over the synodic period 


Therefore, we averaged all parameters obtained from the numerical solution of equation of motion 


(Eq. 59) over the synodic period. The averaged partial derivatives of the disturbing function Rq 


in Eq. (Al) can be calculated numerically. If we know the semimajor axis, eccentricity, longitude 


of pericenter and the resonant angular variable, then we can substitute these values in Eqs. (Al) 


in order to determine the coefficients in Eqs. (16)-(19). The calculation of A 3 , A 2 and Ai using 
Eqs. (16)-(18) gives A* as roots of the cubic equation Eq. (29). If all roots are real, then the 


linearization approximation does not hold well. In the mean motion resonance should be always 
the oscillations of the semimajor axis which maintain the ratio of periods near the ratio of two 
small natural numbers. Therefore, the only usable roots for the librational solution are one real 
number and two complex numbers which are complex conjugate to each other. The imaginary part 
of the two complex A* determines an angular frequency of the libration. 
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6/5, eiim = 0.2472262 


5/4, eiim = 0.2736100 


4/3, eii„ = 0.3107946 





3/2, dim = 0.3690281 


5/3, dim = 0.4140014 


2/1, dim = 0.4811816 





5/2, eiim = 0.5505039 


3/1, dim = 0.5993462 


4/1, dim = 0.6653701 





Fig. 1.— The resonant angular variables obtained from the resonant condition da/dt = 0 for nine 
exact exterior mean motion resonances in the PCR3BP with the PR effect and the radial stellar 
wind. Dust particles in the mean motion resonances with the Earth in a circular orbit around the 
Sun are considered. The resonant angular variables are calculated for the eccentricity of particle 
orbit equal to the universal eccentricity eum. The minimal sizes of the dust particles are determined 
by existences of valid solutions of the linearized system of equation given by Eqs. (12). The star 
in the third plot denotes the position of the evolution depicted in Figure]^ 


It is possible to dehne an exact mean motion resonance. In the exact mean motion resonance 
averaged mean motion of the particle and mean motion of the planet have a ratio that is exactly 
equal to the ratio of two natural numbers. For the exact resonance in the PCR3BP with radiation 
is the semimajor axis equal to a = ap (1-/3)^/^ + Mp)]^^^ [p/iP + Q)]'^^^- Periodicity and 


stability of evolutions in the exact resonances in the PCR3BP with radiation is discussed in Pastor 
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5/2, eiim = 0.5505039 3/1, eii^ = 0.5993462 





4/1, eiim = 0.6653701 




Fig. 2.— The angular frequencies of libration for the dust particles in the mean motion resonances 
with the resonant angular variables shown in Figure The used colors of curves correspond with 
the colors of curves in Figure The star in the third plot denotes the position of evolution depicted 
in Figure In order to show the angular frequencies for the blue curves in Figure the zoom 
squares were used (2/1, 3/1 and 4/1 mean motion resonances). 


(2015). 


An evolution of eccentricity of the particle’s orbit during a capture in an exterior resonance 
approaches the universal eccentricity in the PCR3BP with radiation (Beauge &: Ferraz-Mello||l994 


Liou & Zook 1997). The universal eccentricities are given by 

p + q 


1 - 






= 0 . 


(61) 
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Fig. 3.— A comparison of the results obtained from the numerical solution of the equation of 
motion (solid black curve) with predictions from the analytical theory (dashed red curve) for the 
evolutions of semimajor axis a, eccentricity e, longitude of pericenter oj, resonant angular variable 
a and kh point for a dust particle with /? = 0.025 and = 1 captured in an exterior mean motion 
orbital 4/3 resonance with the Earth under the action of the PR effect and the radial solar wind. 
The initial eccentricity of the particle orbit is equal to the universal eccentricity e\\^ (Eq. 61) and 
as a consequence the kh point librates around an libration center. 


Astrophysical problems with a rotational symmetry around the star lead to the secular time 
derivatives of orbital elements caused by the considered non-gravitational effects that do not depend 
on the longitude of pericenter. The stellar radiation with radial stellar wind also has the rotational 











20 







/ / 

/ / 

/ / 

wr 

\ 

1 - 
1 

\ \ 

\ \ 

\ \ 

■ 1 
/ 

/ 


Fig. 4.— The same parameters in the same PCR3BP with radiation as in Figure are compared. 
However, in this case the dust particle is captured in an exterior mean motion orbital 3/2 resonance. 
The black solid curve corresponds to the evolution obtained from the numerical solution of equation 
of motion. The dashed red curve corresponds to the librational solution given by Eqs. (48). The 
collisions in the kh plane occur on the dashed curves. The initial eccentricity of the particle orbit is 
smaller than the universal eccentricity and the eccentricity asymptotically increases to the universal 
eccentricity. In this case there is not libration of the kh point around an libration center. 


symmetry and {da/dt)EF, in the first equation in Eqs. does not depend on the longitude of 
pericenter w. w has single occurrence in the averaged disturbing function Rq inside the resonant 
angular variable a. Therefore, from a resonant condition da/dt = 0 (substituted in the first equation 
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in 


Eqs. we obtain for the stellar radiation a resonant a which is equal for all values of w. 


Knowing of the resonant a is necessary if we want to determine the angular frequency of 
libration. We calculated the resonant a for nine exact exterior mean motion resonances at the 
universal eccentricities in the PCR3BP with radiation which comprises the Earth and the Sun 
as the two major bodies. The librational solution with complex Aj does not exist for all found 
a. The maximal /? (smallest size of the dust particle) on the red curves in Figure corresponds 
to the maximal /3 for which the librational solution exists on these curves. This rule determines 
the top boundary for (3 in Figure The angular frequency of libration for the problem solved in 
Figure are shown in Figure]^ For /3 = 0 (the PCR3BP) the 2/1, 3/1 and 4/1 exterior mean 
motion resonances in Figure have the asymmetric libration centers (see |Beauge|[l994[ [Beauge fc 


Ferraz-Mello||19M). The asymmetric libration centers are not located at fi = 0 or in the middle of 


the allowed interval for a. In Figure we can see how the red curves rapidly goes to zero at the 
maximal /? after initial increasing at smaller values of (3. The angular frequency for the resonant 
angular variable depicted with the green curves in Figure [^always increases with increasing (3. The 
angular frequency for the blue curves deceases with increasing j3. 


In Figure are compared evolutions of four variables describing the secular orbital evolution of 
the dust particle in the mean motion resonance and evolutions of a kh point for circular Sun-Earth- 
dust system with the solar electromagnetic radiation and the radial solar wind taken into account. 


Coordinates of the kh point are the non-canonical variables k = ecosci and h = esincj (Beauge 


&: Ferraz-Mello 

1993 

Sidlichovsky & Nesvorny 

1994; 

Pastor 

2014 


2014). Collisions of the particle with 


the planet occur on dashed curves in the kh plane. The dashed curves cannot be crossed by the 
kh point during the evolution in a mean motion resonance. Properties of the used dust particle 
are f3 = 0.025 and Q^j. = 1. The particle is placed in an exact exterior mean motion orbital 4/3 
resonance, initially. The initial eccentricity equal to the universal eccentricity enm is used. The 
initial resonant variable is determined from the resonant condition da/dt = 0 (substituted in the 
first equation in Eqs. for the exact resonance with the universal eccentricity. For the numerical 
solution of the equation of motion (Eq. 59) the initial positions of the planet and the particle are 
necessary. The planet was initially located on the positive x-axis and the dust particle was initially 
located in the pericentre. The initial positions of the planet and the particle together with the 
initial resonant angular variable determine the initial longitude of pericenter, according to Eq. ([^, 
as thin = — q/ {p + q) CTm- The black solid curve corresponds to the evolution obtained from the 
numerical solution of the equation of motion. The dashed red curve corresponds to the librational 
solution given by Eqs. (48). Since the longitude of pericentre is advancing we used the term linearly 


proportional to the time (Eqs. [4^ only in the longitude of pericentre (Eqs. 46). As can be seen in 
Figure the librational solution describes the evolution of all parameters sufficiently in this case. 


Figure shows a comparison of the numerical solution of the equation of motion with the 
analytical librational solution for the same dust particle under the action of the same forces as in 
Figure]^ However, the particle is captured in this case in an exterior mean motion orbital 3/2 
resonance. The initial semimajor axis corresponds to the exact resonance. The initial eccentricity 
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of the particle orbit is 0.2. The initial resonant angular variable is 293.8°. The planet was initially 
located on the positive x-axis and the dust particle was initially located in the pericentre. Similarly 
to Figure|^for the initial longitude of pericentre holds Cj^ = — q/ {p + q) o'in- The initial eccentricity 
0.2 for the motion solved in Figure is not equal to the universal eccentricity and the evolution 
of eccentricity monotonically approaches the universal eccentricity, as can be seen in Figure We 
have used the terms varying linearly with time for the eccentricity and the longitude of pericenter 
in this case. Again the librational solution well describes the evolutions of all parameters in Figure 
1^ The position of the kh point decreases more quickly in the analytical solution mainly due to 
differences in the evolutions of the eccentricity between the numerical and the analytical solution. 


6. Discussion 


Figures and depict cases when the linearization approximation holds well. Such cases 
are characterized with small variations in all evolving parameters a, e, oj and a averaged over 
the synodic period. There exist cases when the linearization approximation is not usable. The 
orbital evolutions in mean motion resonances in the PCR3BP with radiation have relatively easy 
characteristic properties. The semimajor axis always oscillates. The eccentricity has oscillations 
with small amplitude that are modulated on the asymptotic approach of the eccentricity to the 
universal eccentricity for the exterior resonances and for the interior resonances the oscillations are 


modulated on a monotonic decrease of the eccentricity to zero ( 

Pastor et al.|2009t 

Pastorp013 

). The 

longitude of pericente usually varies monotonically (but not always, see 

Pastor, Klacka & Komar 


2009). The resonant angular variable has oscillations with an amplitude smaller than vr modulated 


on slower variations. We found that the amplitude of oscillations in the evolution of the resonant 
angular variable has largest influence on the applicability of the librational solution. The librational 
solution is not usable for large amplitude oscillations in the resonant angular variable. We have 
obtained the analytical solution only using an assumption for values of that are unknown. If 
the partial derivatives with respect to time in Eqs. (11) could be exactly determined, then the 


accordance between the analytical and the numerical results should be better. The solution of the 
linearized system of equations (Eqs. can be easily applied for the mean motion resonances 
without influence of the non-gravitational effects (i.e. in the PCR3BP) by setting corresponding 
partial derivatives to zero. 


7. Conclusion 

We have obtained a librational solution for a body with negligible mass captured in a mean 
motion resonance in the PCR3BP with non-gravitational effects included. The solution describes 
the evolutions of over a synodic period averaged semimajor axis, eccentricity, longitude of pericentre 
and resonant angular variable. According to averaged Lagrange’s planetary equations with the 
non-gravitational effects included the evolutions of these four parameters determine the orbital 
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evolution in the mean motion resonance in full. In order to obtain the solution we performed 
a linearization of the averaged equations of motion around an initial point. In the linearization 
partial derivatives of the averaged equations of motion with respect to time are unknown. They can 
be considered analytically as known during the analytical calculation of the librational solution. 
However, for a comparison with real results their values must be found. The calculated analytical 
solution implicates that a linear proportionality to the time can be present in the evolutions of 
the eccentricity, longitude of perihelion and resonant angular variable. A coefficient at the time 
depends on the unknown partial derivatives with respect to time. We assumed that the value of 
the coefficient can be approximated by over a libration period averaged slope of the eccentricity, 
longitude of pericenter and resonant angular variable. The libration period can be determined 
without knowing of the partial derivatives with respect to time. We compared the analytical and the 
numerical results for the PR effect and the radial stellar wind used as the non-gravitational effects. 
Agreement between the compared results was good. The librational solution can be used in the 
cases when the librational amplitude is small i.e. in the cases when the linearization approximation 
can be used. The comparson of analytical and numerical results for very small libration amplitudes 
is difficult because the small libration can be lost in the numerical averaging process. 


A. Constant coefficients 


This appendix presents coefficients for linearized system of equations (Eqs. 12) describing the 
orbital evolutions of the dust particles captured in the mean motion resonances in the PCR3BP 
with radiation. 
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